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ABSTRACT 

We study the evolution of globular clusters with mass spectra under the influence of 
the steady Galactic tidal field, including the effects of velocity anisotropy. Similar to 
single-mass models, velocity anisotropy develops as the cluster evolves, but the degree 
of anisotropy is much smaller than isolated clusters. Except for very early epochs of 
the cluster evolution, nearly all mass components become tangentially anisotropic at 
the outer parts. We have compared our results with multi-mass, King-Michie models. 
The isotropic King model better fits to the Fokker-Planck results because of tangential 
anisotropy. However, it is almost impossible to fit the computed density profiles to the 
multi-mass King models for all mass components. Thus if one attempts to derive global 
mass function based on the observed mass function in limited radial range using multi- 
mass King models, one may get somewhat erratic results, especially for low mass stars. 
We have examined how the mass function changes in time. Specifically, we find that 
the power-law index of the mass function decreases monotonically with the total mass 
of the cluster. This appears to be consistent with the behaviour of the observed slopes 
of mass functions for a limited number of clusters, although it is premature to compare 
quantitatively because there are other mechanisms in contributing the evaporation of 
stars from the clusters. The projected velocity profiles for anisotropic models with the 
apocenter criterion for evaporation show significant flattening toward the tidal radius 
compared to isotropic model or anisotropic model with the energy criterion. Such a 
behaviour of velocity profile appears to be consistent with the observed profiles of 
collapsed cluster Ml 5. 
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1 INTRODUCTION 

The study of evolution of globular clusters has a long his- 
tory. It would be most desirable to use direct A^-body calcu- 
lations in order to follow the evolution of realistic globular 
clusters realistically, but we do not have enough computing 
power to do that for A*' comparable to the number of stars in 
real globular clusters at this moment. Therefore we have to 
rely on approximate techniques to understand the dynam- 
ics of star clusters. The Fokker-Planck equation has been 
widely used for the study of globular clusters, but it requires 
many simplifying assumptions. The assumption of velocity 
isotropy is usually employed because it allows fast and accu- 
rate integration of the Fokker-Planck equation (Cohn 1980) . 
However, high degree of velocity anisotropy is known to be 
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generated by two-body relaxation in the outer parts of iso- 
lated clusters (e.g., Spitzer 1987). The observational prop- 
erties of the clusters depend on the degree of anisotropy 
and correct treatment of anisotropy becomes an important 
issue. The development of velocity anisotropy was investi- 
gated in some recent works which used various numerical 
methods: e.g., direct A'^-body calculations by Giersz & Heg- 
gie (1994), gaseous model calculations by Louis & Spurzem 
(1991) and Spurzem (1996), and Fokker-Planck calculations 
by Cohn (1979), Drukier et al. (1999) and Takahashi (1995, 
1996, 1997). Most of these calculations confirmed the earlier 
expectations of generation of radial anisotropy in the outer 
parts of the clusters, although detailed behaviour depends 
somewhat on the numerical methods. 

Actual clusters are embedded in the Galactic tidal field 
which should impose a finite boundary to the cluster. The 
introduction of tidal boundaries should have significant ef- 
fects on the velocity anisotropy. Recently Takahashi, Lee 
& Inagaki (1997; abbreviated as TLI hereafter) performed 
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anisotropic Fokker-Planck simulations for single mass clus- 
ters and demonstrated that radial anisotropy grows up to 
the time of core collapse, but the degree of anisotropy be- 
gins to decrease to smaller values as the cluster expands and 
loses significant amount of mass. They showed that there is 
even tangential anisotropy very close to the tidal boundary 
in the late phase of the cluster's evolution. The density pro- 
files of the clusters can be reasonably fitted by isotropic King 
models for most of the time except for the brief phase around 
the core collapse, when the degree of radial anisotropy nearly 
reaches its maximum. 

The behaviour of the anisotropy for isolated clusters 
with mass spectra has some interesting features. Takahashi 
(1997) showed that the degree of anisotropy for different 
mass components behaves differently. There is always ra- 
dial anisotropy for low mass components, but tangential 
anisotropy can develop for high mass components at early 
evolutionary stages when the temperature of heavy stars de- 
creases rapidly toward the equipartition. During this 'cool- 
ing' phase, the radial velocity dispersion decreases more 
rapidly than the tangential one. This leads to the devel- 
opment of tangential anisotropy among high mass stars. 

As in the case of single mass models, extension of these 
studies into tidally limited clusters should provide useful in- 
sight on the evolution of realistic clusters. Based on the expe- 
rience of single mass models, one obviously expects that the 
velocity anisotropy will be greatly affected by the presence 
of the tidal field. Some of the essential points have been clar- 
ified by the analyses of TV-body results (Vesperini & Heggie 

1997) , but one should remember that the A''-body calcula- 
tions are done for small N cases. Many of the behaviour of 
self-gravitating systems does not allow simple scaling (e.g., 
Portegies Zwart et al. 1998, Takahashi & Portegies Zwart 

1998) . We make detailed calculations for the evolution of 
globular clusters with these effects in order to make quanti- 
tative assessments for the long term evolution of star clusters 
with the tidal cut-off and mass spectrum. 

The study of the evolution of tidally limited clusters 
has important areas of application. The mass function of 
globular clusters is difficult to obtain because of large dy- 
namic range in brightness of individual stars as well as large 
variation of stellar densities. It is very difficult to determine 
the global mass function observationally, but there have been 
several studies to determine the global mass function based 
on locally measured mass functions in limited areas of globu- 
lar clusters (Richer et al. 1990, 1991 and references therein). 
It is absolute necessity to obtain the global mass function 
if one wants to study factors determining the mass func- 
tion of clusters. The conversion of the local mass function 
to the global one is not a trivial task. The general prac- 
tice employed widely is to fit the observation to multi-mass 
King models (e.g., Meylan 1988) or to assume that the mass 
function in low mass stars in the outer parts are close to the 
global one (Richer et al. 1991). However, there is no guar- 
antee that these processes give a right answer. The detailed 
studies of evolving clusters with the tidal boundary can pro- 
vide an important check of the above mentioned practices. 

The mass function changes with time because the rate 
of escape of stars depends on the mass of the star. The high 
mass stars tend to spiral into the central parts by dynam- 
ical friction, and consequently low mass stars are preferen- 
tially removed from the tidally limited clusters. The slope 



of mass function thus becomes less steep as the cluster loses 
mass. Even if one gets the global mass function through care- 
ful observations and correct conversion process, the present 
day mass function (PDMF) may be significantly different 
from the initial mass function (IMF). The detailed studies 
of evolving star clusters could provide an important frame- 
work of estimating the relation between PDMF and IMF. 

In the present paper, we make detailed studies of the 
the evolution of tidally limited multi-mass clusters including 
the effects of velocity anisotropy. We pay special attention to 
the behaviour of mass function during the course of dynam- 
ical evolution. We examine the adequacy of the process of 
converting the local mass function to the global mass func- 
tion using multi-mass King model fitting. We then analyze 
how the mass function changes with time and location. 

This paper is organized as follows. In §|^, we describe 
our models for the evolution of globular clusters. In §|^, we 
introduce the idea of half-mass concentration and discuss its 
relation to the development of anisotropy and to mass loss 
rates. InjM we describe initial conditions of multi-mass mod- 
els. In §b[we present the results of the model calculations. 
The relationship between local and global mass function and 
the temporal evolution of mass function are examined in §^ 
We discuss the velocity dispersion profiles in § |^. The final 
section summarizes basic results. 



2 FOKKER-PLANCK MODELS 

We use the anisotropic Fokker-Planck code developed by 
Takahashi (1995). It is based on the direct integration 
method of the Fokker-Planck equation in (E, J) space (Cohn 
1979), where E = + (j){r) is the energy of a star per 
unit mass and J = rvt is the angular momentum per unit 
mass. The extension to multi-mass clusters is described in 
Takahashi (1997) and the treatment of tidal boundary can 
be found from TLI. TLI used two different criteria for the 
escape of a star from the cluster: apocenter criterion and en- 
ergy criterion. The apocenter criterion assumes that a star 
escapes if the apocenter distance of a star with {E, J) be- 
comes greater than tidal radius rj. This criterion preferen- 
tially removes the stars with low J and thus suppresses the 
development of radial anisotropy. The energy criterion as- 
sumes that a star escapes if E > Et where Et — (l>{rt)- 
The radial anisotropy is also suppressed in this criterion and 
tangential anisotropy can develop (see TLI). The general be- 
haviour of the cluster evolution is not much dependent on 
the choice of escape criteria for the models of TLI, but de- 
tails of the clusters vary with the treatment of stellar evapo- 
ration. Takahashi & Portegies Zwart (1998) have shown that 
the apocenter criterion gives better results when compared 
with realistic A''-body calculations. In the present study, we 
examine the results of our calculations with both of the cri- 
teria, but more emphasis is given to the calculations with 
the apocenter criterion. 

Hereafter, we denote anisotropic models with the apoc- 
enter criterion as Aa models and anisotropic models with 
the energy criterion as Ac models. The capital A stands 
for anisotropic models, and the lower cases a and e stand 
for apocenter and energy criteria, respectively. Similarly 
isotropic models are denoted as le models, but la models 
do not exist. 
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3 HALF-MASS CONCENTRATION, 

VELOCITY ANISOTROPY, AND THE MASS 
LOSS RATE 

Before going into descriptions of the results of multi-mass 
models, we examine the effects of the degree of the concen- 
tration of a cluster on the development of velocity anisotropy 
and on the cluster mass loss rate. We use single-mass models 
in this section because these models are simple, but the con- 
clusions of this section are also true for multi-mass models, 
as we show later. 

TLI used the initial model of a single-mass Plummer 
model with the cutoff radius of 32 in Plummer model units. 
The results of their simulations show that the Aa and Ae 
models lose mass faster than the le model (see Fig. 2a of 
TLI, though the Ae model is not shown in it). The evapora- 
tion time of the le model is about 12% longer than that of 
the Aa model, and the latter is about 10% longer than that 
of the Ae model. The biggest difference in the mass loss rate 
between the anisotropic models and the isotropic model ap- 
pears in the pre-coUapse and the early post-collapse phases. 
The mass loss rates of the different models are similar in the 
late post-collapse phase. 

In §1^ we show the results of simulations which start 
with multi-mass King models. As shown in Fig. ^ the le 
model and the Ae model are very similar in the evolution of 
the total mass. The Aa model loses mass slower than the le 
model. This trend is opposite to that seen in the simulations 
of TLI. While Fig. § shows the results for a. Wo = 3 King 
model, the same trend is seen for King models with other 
Wq. Why does the opposite trend appear? The change from 
the single-mass model to the multi-mass model is not the 
reason of the change of the trend. Rather, it is found that 
the trend changes with the concentration of a cluster to the 
tidal cutoff radius. 

First we discuss what parameter is useful in measur- 
ing the degree of the concentration of a cluster. For King 
models, the concentration is usually measured in terms of 
the concentration parameter c = logjQ(rt/ro), where vq is 
the King core radius (King 1966). The parameter c is useful 
in measuring the concentration of the core. However, when 
discussing the variation of the total mass, c is not always 
an appropriate parameter, because the fraction of the core 
mass is very small in large c clusters. It may be more ap- 
propriate to use the half-mass radius rh instead of the core 
radius vq. Thus we define half-mass concentration 



Ch = n/^h, 



(1) 



(note the logarithm is not taken here). 

King models form a one-parameter family parameter- 
ized by the concentration c or equivalently by the dimension- 
less central potential depth Wo; c monotonically increases 
with Wo (King 1966). On the other hand, Ch does vary mono- 
tonically with Wo- The relation among Wo, c, and Ch is given 
in Fig. hi for 1 < Wo < 14. In contrast to that the core con- 
centration c varies by nearly three orders of magnitude, the 
variation of the half- mass concentration Ch is within a factor 
of three for this range of Wo - The half-mass concentration Ch 
increases with Wo for Wo = 1 to 8, then it decreases until Wo 
reaches and increases again. For the range of 1 < Wo < 14, 
Ch reaches the maximum of about 10 at Wo — 8. The initial 
model of TLI has Ch = 25, therefore it is far more con- 
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Figure 1. The behaviour of Ch, c, and ^as for King models 
against Wo- Note that cj, has a maximum, and ^as has a broad 
minimum near Wo — 8. 



centrated than any King models in terms of the half-mass 
concentration. 

In order to investigate the effects of the half-mass con- 
centration on the cluster evolution, we carry out simulations 
which start with single-mass Plummer models with different 
cutoff radii which correspond to Ch =10, 20, and 30. In these 
simulations the number of particle is set to 10*, which is 
specified to include the effects of three-body binary heating. 
Fig. ^ shows the evolution of the total mass for the Ae and 
le models. For the large Ch initial models, the Ae model loses 
mass much faster than the le model at the early evolution- 
ary phase as in TLI's simulations. (The core collapse occurs 
at t =16-18 trh,o in all the models.) As Ch decreases, the dif- 
ference between the Ae and le models decreases and finally 
almost diminishes for Ch = 10. Fig. ^ shows the evolution 
of the anisotropy /3 at the 90%-mass radius. The anisotropy 
parameter l3 is defined as follows; 



(2) 



where ot and ar are one-dimensional tangential and ra- 
dial velocity dispersions, respectively. The development of 
the radial anisotropy during the pre-coUapse phase becomes 
strongly suppressed as Ch decreases. This is simply because 
radial-orbit stars can easily escape from a small tidal radius 
cluster. As TLI discussed, the emergence of a large num- 
ber of radial-orbit stars is the reason for faster mass loss in 
the anisotropic models. For low Ch clusters, since there is 
almost no room for the radial anisotropy to develop, the Ae 
model and the le model do not differ very much in the evo- 
lution of the total mass. Since King models have such low 
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Figure 2. The evolution of mass as a function of time for Ae 
models (dotted lines) and le models (dashed lines). The initial 
models are single-mass Plummer models with three different tidal 
cutoff radii: Ch = rt/rh =10, 20, 30, from left to right. As Ch 
decreases the mass loss becomes faster and the difference between 
the Ae model and the le model becomes smaller. The core collapse 
occurs at t =16— IStrh o iu all the models. 



Ch (ch < 10), the difFerence among the three models shown 
in Fig. ^ can be understood from this viewpoint. The mass 
loss of the Aa model should be always slower than that of 
the Ae model because of the apocenter criterion (see TLI). 

Note that we have assumed that the initial tidal radius 
is equal to the limiting radius of King models. This is a con- 
ventional way for setting up initial conditions of simulations, 
but not the only one; these two radii can be different. 

Finally in this section, we comment on the relation be- 
tween the mass loss rate and the half-mass concentration 
for King models. Johnstone (1993, figure 2) and Gnedin & 
Ostriker (1997, figure 6) found that for King models the 
mass loss rate per half-mass relaxation time (exactly speak- 
ing, the inverse of the evaporation time in units of the ini- 
tial half-mass relaxation time in Gnedin & Ostriker) has a 
minimum around c =1.5-2. Johnstone (1993) calculated the 
evaporation rate for King models using the local Fokker- 
Planck equation, and Gnedin & Ostriker (1997) performed 
isotropic orbit-averaged Fokker-Planck simulations with the 
initial conditions of King models. Their findings are quali- 
tatively consistent with a simple estimation of the mass loss 
rate based on only the half-mass concentration c^. Naively 
one may think that the mass loss rate decreases as Ch in- 
creases. In fact Ch has the maximum at Wo ~ 8 (c ~ 1.8), 
where Johnstone (1993) and Gnedin & Ostriker (1997) found 
the minimum evaporation rate. In addition Johnstone (1993) 
found a local maximum of the mass loss rate at c ~ 2.6. This 
may correspond to the fact that Ch has a local minimum at 
c ~ 2.7 {Wo — 12) as shown in Fig. ^. As an additional 
reference, we show the mass loss rate per half-mass relax- 
ation time ^AS calculated from the modified Ambartsumian- 
Spitzer formula (Spitzer 1987, p. 57; see also TLI) for King 
models in Fig. |l|. This estimation is even quantitatively not 
very different from the result of more detailed calculations 
for W^o < 8 by Johnstone (1993). 
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Figure 3. The evolution of the anisotrony /3 at the 90%-mass 
radius for the Ae models plotted in Fig. ^ The initial develop- 
ment of the radial anisotropy becomes strongly suppressed as cji 
decreases from 30 to 10. 



4 MULTI-MASS MODELS 

There are growing evidences that the mass function in 
clusters as well as in the field is not a simple power-law, 
but should be approximated as composite power-laws (e.g., 
Kroupa, Tout & Gilmore 1993). However, we use simple 
power-laws in the present study of general evolution of 
clusters, for simplicity. The adoption of multiple compos- 
ite power-laws is trivial task, but requires specification of 
more model parameters. The assumed initial mass function 
(IMF) has the following form: 



N(m)dm = Cm "dm, 



<m< rrin 



(3) 



where C and a are constants. As in Takahashi (1997), we use 
the discrete mass components. The i-th mass bin for equal 
A log m interval can be obtained by 

m, = m^in(^^) , (i = l,2,...,if), (4) 

V mmin / 

where K is the number of components. The parameters de- 
termining the mass function are a and mmax/T7^min). 
We use a few values of a and ^ for our model calculations. 

The issue of the effects of stellar evolution naturally 
comes out for multi-mass models if one extends the mass 
function to high mass stars. Obviously the stellar evolution 
should have been the most important process in determin- 
ing the early phase of the dynamical evolution of globular 
clusters. However, we will ignore the stellar evolution in this 
paper because it adds one more complexity to the models. 

Without stellar evolution, the clusters usually undergo 
core-collapse before losing significant amount of stars unless 
the initial cluster has very low central concentration, because 
the time scale for the core collapse is generally much shorter 
than the evaporation time. Note that the core collapse time 
is of order of lOtrh.o while the evaporation time is typically of 
order of 40trh,o, where frh.o is the initial half-mass relaxation 
time for single mass models. The core collapse is very much 
accelerated for multi-mass clusters, and the discrepancy be- 
tween the evaporation time and the core collapse time could 
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become even larger. Therefore one needs to supply the mech- 
anism that allows the evolution of cluster beyond the core 
collapse. We assumed that the heating is provided by bina- 
ries formed by three-body processes because it is easiest to 
implement in the framework of Fokker-Planck method. The 
evolution of parameters for the central parts such as den- 
sity and velocity dispersion should depend on the detailed 
mechanism for driving the post core-collapse evolution, but 
most of the global properties do not sensitively depend on 
the heating mechanism. For our purposes, three-body binary 
heating is sufficiently simple and easy to be used as heating 
mechanism. The heating formula for multi-mass clusters is 
taken from Lee, Fahlman & Richer (1991). 



5 RESULTS OF MULTI-MASS MODEL 
CALCULATIONS 

5.1 Initial Model Parameters 

We have used various initial models. The characteristics of 
the initial models are summarized in Table 1. In all mod- 
els, we have assumed that the shape of density profiles of 
different mass components are the same with the different 
amplitudes dictated by the assumed mass functions. The 
velocity dispersions are assumed to be the same for all com- 
ponents, are assumed to have the same density and velocity 
profiles. 

The behaviour of the cluster up to the core collapse does 
not depend on the initial total mass or size of the cluster. 
But the amount of binary heating depends on the number 
of stars in the cluster. We have assumed that the mass of 
the cluster is 6 x 10* Mq for the model R, and we have 
fixed the total number of stars at 5 x 10* for other models. 
But this number only controls the central density during 
the postcoUapse phase. Most of our results does not really 
depend on the mass of the cluster, except for the surface 
density profiles. If we used higher initial mass, the core radii 
at a given epoch (measured by the M/Mq) would be smaller. 
The minimum mass rrimin of the mass function was assumed 
to be 0.1 M0 in all models. 

Throughout this paper, we use r to represent the dis- 
tance from the center of the cluster, and R to represent the 
projected distance. 

5.2 A Reference Model 

We first describe the general behaviour of Reference Model 
(Model R) whose initial parameters are given in Table 1. 
This set of parameters is used by Heggie et al. (1998) to 
compare the outcomes of various computational methods. 
We will consider this initial condition as a reference model. 
The number of mass component was 8. 

5.2.1 Global behaviour 

The evolution of the total mass is shown in Fig. ^ with dif- 
ferent criteria for ejection of stars. Unlike single component 
models, the mass loss rate is not a constant, but increases 
with time. The increase of M is due to the fact that the mass 
function becomes flatter as the cluster loses mass. (The mass 




Figure 4. The evolution of mass of the model R of three different 
assumptions (Aa, Ae, and le). The mass loss rate is nearly a 
constant for single component models, but increases with time 
in our multi-mass models. This is due to the fact that the mass 
function becomes flatter as the cluster loses mass. The flatter 
mass function means higher mass loss rate. 

loss rate becomes larger for clusters with flatter mass func- 
tion.) As discussed in §3, the apocenter criterion (Aa) gives 
longer lifetime than the energy criterion (Ae), but it is not 
a significant effect. For comparison, we have also shown the 
mass evolution for the isotropic model (Je) with the same 
initial conditions. The time of maximum collapse is shown 
as dots in this figure. 

5.2.2 Degree of anisotropy 

The behaviours of the degree of velocity anisotropy l3 for 
the apocenter criterion are shown in Figure The upper, 
middle and lower panels show P{r) at the time M/Mo=0.7, 
0.5 and 0.2, respectively. The same figure as Fig. ^ for the 
energy criterion is shown in Fig. ^ 

Note that the sign of (3 can be both positive (radial 
anisotropy) and negative (tangential anisotropy). If a clus- 
ter is isolated, strong radial anisotropy develops in the halo 
for all components, as shown by Takahashi (1997). However 
the present model cluster is in a rather strong tidal field: 
the initial King model with Wo ~ 3 has Ch = 3.7. Therefore, 
from the discussion in §^ we may expect that the develop- 
ment of such radial anisotropy is strongly suppressed in the 
present case. From the profiles shown in Fig. |^ {Aa model), 
in fact, only weak radial anisotropy is seen for the high mass 
components at intermediate radii. The low mass components 
tend to have tangential anisotropy at all radii. A similar 
trend is seen in Fig. ^ {Ae model). This is in sharp contrast 
with isolated clusters where strong radial anisotropy devel- 
ops among low mass components (Takahashi 1997). In the 
top panel of Fig. 0, the high mass components also have tan- 
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Table 1. Initial Parameters for Multi-Mass Models 



Model 


Profile 


a 




Comments 


R 


King, Wo = 3 


2.35 


15 


reference 


Kl 


King, Wo = 3 


2.35 
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low fl 


K2 


King, Wo = 8 


2.35 


15 


high Wo (ch) 
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Plummer, cj^ = 20 
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Figure 5. The dependence of /3 on r at three different epochs: 
M/Mo=0.7, 0.5 and 0.2 for R models with the apocenter condi- 
tion. The core-collapse occurs when M/Mo ~ 0.70. Note that the 
low mass components show more tangential anisotropy than high 
mass stars. This is opposite to the isolated cluster models. 
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Figure 6. The dependence of /3 on r at three different epochs: 
M/Mo=0.7, 0.5 and 0.2 for R models with energy condition. Note 
that the core-collapse occurs when M/Mo ~ 0.57. The velocity 
distribution is forced to become isotropic at tidal boundary in the 
energy condition. 



gential anisotropy. This is due to the initial cooling of the 
high mass components (Takahashi 1997). The velocity dis- 
tribution becomes mostly tangential near the tidal boundary 
regardless of the mass. The tangential anisotropy was also 
noted by Oh & Lin (1992), and Kim & Oh (1999) from their 
hybrid calculation of N-body and Fokker-Planck. 

The behaviour of /3(r) becomes different for apocenter 
and energy criteria in the outer parts, because the energy 
criterion forces /3 at tidal boundary to become zero (i.e., 
isotropic) . 

5.2.3 Density profiles 

The profiles of surface densities for different mass compo- 
nents are shown in Fig. |^ for a model with the apocenter 
condition at the epochs of M/Mo = 0.7, 0.5 and 0.2. Fig. | 
shows the density profiles for a model with the energy condi- 
tion. The density profiles at M/Mo=0.5 and 0.2 are similar 
for the apocenter condition and the energy condition. The 



big difference at M/Mo— 0.7 is purely due to the fact that 
the core-collapse occurs when M/Mo ~ 0.7 for the Aa model 
while it occurs somewhat later for the Ae model. Therefore, 
the difference in the degree of anisotropy at outer radii does 
not play important role in determining the density distribu- 
tion. 

The profiles at M/A/o=0.5 and 0.8 are compared with 
the isotropic model in Fig. ^. The radial profiles of the to- 
tal density of isotropic and anisotropic models have very 
little difference, but the detailed distribution of individual 
mass components could be significantly different for these 
two models. From these comparisons, we expect that there 
will be some difference in the behaviour of the mass function 
in the lower mass parts at late phase of the evolution. 

In Fig. |l^, we tried to compare the surface density 
profiles of the R/Aa model at M/Mo=0.5 with the best- 
fitting multi-mass King models having the same mass func- 
tion. Multi-mass King models have the following distribution 
function (e.g., Binney & Tremaine 1987): 
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Figure 7. Surface density distribution of individual mass components at three epochs for the R model with apocenter criterion (Aa). 
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Figure 8. Surface density distribution of individual mass components at three epochs for the R model with energy criterion (Ae). 
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Figure 9. Comparison of density profiles of isotropic (/e: thin lines) and anisotropic (Aa: thick lines) Fokker-Planck results at M/Mo=0.5 
and 0.2. 



ME) = a [eM'E/af) - exp{-Et/a^ 



a = 1,2, 



where at satisfies the relation rriiaf = rnja'j for all i and j. 
The constants d should be determined by the given mass 
function. As in Fig. ^, the profiles of density of low mass 
stars do not fit King models: if the distribution of combined 
density is to be reproduced, the central densities of low mass 
components of King model are higher than that generated 
by Fokker-Planck results. This means that the lower mass 
stars are more distributed at the outer parts than expected 
from King models. The mass segregation effect is greater 
than that predicted by multi-mass King models. This is the 
same trend that was seen from the comparisons between 
anisotropic and isotropic models. 

For the purpose of application to more realistic cases, 
we also attempted to fit the surface density profiles of vis- 
ible components only, assuming that the highest mass bin 
(m = 1.3 Mq) represents the remnant stars such as neutron 
stars or massive white dwarfs. This is still far from more re- 
alistic situation because the extension of simple power-law 
to massive remnant will give results in too much contribu- 
tion from these stars. Nonetheless, such a comparison could 
be worthwhile in understanding the observed data. The re- 
sult is shown in Fig. The fitting model is the same as 
shown in Fig. |l^, but the difference between the King model 
and the Fokker-Planck model is more pronounced. Choosing 
other parameters for King models would not help to improve 
the fitting: the difference simply is due to the fact that the 
distribution of stars are substantially different between King 
models and Fokker-Planck results. Because of stronger mass 
segregation in Fokker-Planck results, the fitting of profiles in 
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Figure 10. Comparison of the surface density profiles of 
Fokker-Planck model at M/Mo=0.5 and the best-fitting 
mass King model. 
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Visible only, M/M„ = 0^5 
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Figure 11. Same as Fig. |lO|, but we have attempted to fit only the 
visible components, assuming that the highest mass component 
is compact stars that have no luminosity. 



the intermediate region gives significant differences in outer 
parts. 

The fitting does not improve when one attempts to use 
anisotropic King models (King-Michie models), because the 
King-Michie models always have radial anisotropy regard- 
less of mass, while the Fokker-Planck results give both ra- 
dial and tangential anisotropy among different mass bins. 
The isotropic King models are usually better in reproducing 
the Fokker-Planck results. This is somewhat contradictory 
with the single component models: during the early phase 
of the evolution, King-Michie models better represent the 
anisotropic Fokker-Planck results of TLI than isotropic King 
models. 



5.3 High Ch cases 

As described in §^for single-mass models, the development 
of radial velocity anisotropy in the outer halo is strongly 
suppressed in low half-mass concentration (ch) clusters. To 
confirm this for multi-mass clusters we now examine the 
results for different initial Ch. 

Fig. |l^ shows the radial profiles of anisotropy /3 at se- 
lected epochs for a reference initial condition except for 
Wo = 8. Initially this model has a larger cii(= 9.6) than 
the model discussed in § (ch — 3.7). Comparing Figs. ^ 
and |l^, however, we find that the behaviours of anisotropy 
are not very different between the two models. For exam- 
ple, the maximum value of f} is about 0.3 at the epoch of 
M/Mo = 0.5 for both models. This is not surprising because 
we may well expect that Ch for the King model of Wo = 8 is 
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Figure 12. Same as Fig. Q but for a model for the initial con- 
ditions of the King model with Wo = 8 (i.e., Kl which have 
Cjj = 9.6), a = 2.35, and /i = 15. The core-collapse occurs when 
Af/Mo=0.99. 



still not large enough for strong radial anisotropy to develop 
(ch < 10), if we remember the results presented in § |^. 

Fig. |l^ shows P profiles for an Aa model for the initial 
conditions of a Plummer model with Ch = 20, a = 2.35, 
and fi = 15 (model PI). In this case, as expected, stronger 
radial anisotropy (maximum /3 ~ 0.5) appears and is kept 
for a longer time than in the previous two cases. 



6 THE MASS FUNCTION 

It is observationally challenging task to obtain the mass 
function of globular clusters. Globular clusters are known as 
oldest Galactic objects and should possess very important 
memory of the early phase of the universe. For example, 
the IMF of globular clusters could provide a glimpse of the 
star formation process during the phase of very low heavy 
elements. 



6.1 Evolution of The Mass Function 

The loss of mass occurs near the tidal boundary. The high 
mass stars gradually spiral into the central parts via dy- 
namical friction and low mass stars are preferentially re- 
moved from the cluster. Therefore the shape of the mass 
function changes with time. This phenomenon has been ob- 
served from many numerical simulations (e. g., Chernoff & 
Weinberg 1990; Lee, Fahlman & Richer 1991; Lee & Good- 
man 1995; Vesperini & Heggie 1997). The evolution of the 
mass function at several different epochs are shown in Fig. 
[l4| . The epochs are chosen according to the mass of the clus- 
ter: M{t)/Mo=Q.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2 and 0.1. This 
figure clearly shows the flattening of the mass function as 
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Figure 13. Same as Fig. ^ but for a model for the initial condi- 
tions of the Plummer model (model PI) with = 20, a = 2.35, 
and fi = 15. The core-collapse occurs when M/Mo=0.99. 



the cluster loses mass. Also shown in this figure are the mass 
functions within the half-mass radius as dotted lines. 

The shape of the mass function deviates from the simple 
power-laws as a result of the selective mass loss: the aver- 
age slope of the mass function decreases with time, and the 
slope of the mass function is steeper at lower mass than at 
higher mass. Although the mass loss rate is larger for the 
lower mass, it is not sensitive enough to flatten toward the 
lower mass. Such a behaviour (i.e., steeper slope at lower 
mass parts) appears to be inconsistent with the observa- 
tional data for a number of globular clusters (e.g., Chabrier 
& Mera 1997). The initial mass function of the globular clus- 
ters should not be simple power-laws since the dynamical 
evolution only makes the situation worse, but it is likely to 
flatyten (or even decreasing) toward the lower mass. 

The evolution of mass function as a function of time 
can be measured by the variation of the power-law index. 
Suppose that the mass function can be approximated as a 
power- law with index a as shown in equation (^). As seen 
from Fig. it is clear that the mass function does not 
follow a single power law as the cluster loses mass. Thus 
we have computed instantaneous a at two different masses: 
m = 0.2 Mq and 0.55 Mq. The behaviour of these indices is 
shown in Fig |l| as a function of M{t) /M (0). The solid lines 
represent a(0.2 Mq) and the dotted lines are for a(0.55 Mq). 
Since the mass function depends on the location as well, we 
have computed q's within the entire cluster (thick lines) and 
within the half-mass radius (thin lines). 

The trends toward the flatter mass function with time 
are clearly shown in this figure. The global mass function 
changes slowly in the early phase of the evolution, and rather 
abruptly in the late phase. The mass function within the 
half-mass radius changes more rapidly than the global mass 




1 F 



0.1 =- 



0.01 ^ 



0.001 



0.0001 



Figure 14. The mass function for R/Aa model at epochs of 
M/Mo=0.8, 0.7, 0.6, 0,5 0.4, 0.3, 0.2 and 0.1 (from steep to flat 
lines). The solid lines represent the global mass function and the 
dotted lines are for the mass functions within the half-mass ra- 
dius. 



function. The model with smaller fi gives more rapid vari- 
ation in a with time than that with higher fi, as shown in 
Figs. ^ and |l|. We also note here that the behaviour of a 
for the isotropic Fokker-Planck model is very close to that 
of the anisotropic model. 

The variation of mass function observed by Richer et al. 
(1991) could be most attributed by the dynamical evolution 
since the mass function tends to be flatter for clusters with 
short relaxation time. Clusters with short relaxation time 
are likely to have lost significant fraction of the initial mass, 
mostly in the form of low mass stars. Richer et al (1991) 
claim that there exists a tight relation between the disrup- 
tion time Td and the power-law index of the mass function. 
Although the amount of uncertainties is large, the trend of 
flattening of mass function for clusters with significant mass 
loss appears to be real. 

6.2 The Relation between Local and Global MF 

The accurate photometry over any entire globular cluster is 
very difficult. The luminosity function (and thus mass func- 
tion) is usually measured in limited range of the cluster. 
Since the mass function varies with the location if the clus- 
ter has undergone significant evolution, it is important to 
understand the relationship between the locally determined 
mass function and the 'global' mass function. 

We have shown the variation of a as a function of pro- 
jected radius for three different epochs in Fig. ^ for 0.2 M© 
(dotted lines) and 0.55 Mq (solid lines). Also shown as hor- 
izontal lines are the power-law indices of the global mass 
function at that epoch. We have seen from Fig. uA that the 
global mass functions measured at 0.2 and 0.55 Mq are 
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Figure 15. The evolution of power-law indices of the global (thick 
lines) and half-mass radius mass function for R/Aa model, as a 
function of cluster mass. The indices are computed for low mass 
(m=0.2 Mq: solid lines) and intermediate mass (m=0.55 Mq : dot- 
ted lines). The mass function is less steep at higher mass than at 
lower mass. 
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Figure 17. The variation of power-indices of mass functions over 
the projected radius for three different epochs: M=0.9, 0.51 and 
0.11 Mo for /i = 7 model with Salpeter IMF (i.e., model Kl/Aa). 
The solid and dotted lines represent for m=0.2 and 0.55 Mq, 
respectively. The horizontal lines are the indices for the global 
mass function. 
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Figure 16. Same as Fig. |l5| except for IJ, = 7 (i.e., model Kl/Aa). 
Note that the difference between power indices at 0.2 and 0.55 Mq 
are small compared to the fi = 15 case. 



rather similar. However, the local mass function changes by 
a great amount. The 0.2 Mq mass function varies slowly over 
radius, but 0.55 Mq one varies a lot. This is due to the fact 
that the high mass stars are concentrated toward the cen- 
tral parts. Thus one has to be very careful in determining 
the mass function from the observations. The power-law in- 
dex of the global mass function at 0.2 Mq roughly coincides 
with the local value near the half-mass radius. The mass 
function near the half-mass radius appears to represent the 
global mass function, when measured at the lower mass. This 
is also seen in N-body calculations by Vesperini & Heggie 
(1997). The power-law index at rn for 0.55 Mq, however, 
deviates significantly from the global value. Therefore, the 
determination of global mass function remains observation- 
ally challenging task. 



7 VELOCITY PROFILES 

So far we have examined density profiles and mass functions. 
We should be able to obtain more detailed dynamical infor- 
mation by measuring the velocity dispersion profiles of the 
stars in the cluster. In principle, both projected and trans- 
verse velocities can be measured if the astrometric accu- 
racy becomes of order of 10~^ arcseconds, as the astrometric 
project of GAIA is targeting. However, only the projected 
velocity information is available for most of the clusters in 
the Galaxy. We now discuss the prediction of projected ve- 
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Figure 18. Projected velocity dispersion profiles for R/Aa and R/Ie models at three different epochs: M/Mq = 0.7, 0.5, and 0.2. The 
post collapse expansion begins at near M/Mq « 0.7. Note that the velocity profiles are much flatter for the anisotropic model with 
apocenter criterion than the isotropic model. 



locity dispersions of the cluster stars by our Fokker-Planck 
models. 

In Fig. |l^, we have plotted the projected velocity dis- 
persion profiles for the model R/Aa at Af/Mo=0.7, 0.5 and 
0.2 for four different mass components, together with the 
corresponding le at the same epoch (measured by M/Mo). 
The maximum collapse occurs near M/Mq— 0.7. The unit of 
velocity in this plot is GMq/tq. The behaviour of veloc- 
ity profiles for R/Ae model was not shown here, but we note 
that the velocity profiles of such a model behave similarly to 
those of isotropic model. This may be due to the similarity 
in tidal boundary condition. 

The Aa model predicts the velocity profiles flatter than 
the isotropic model. Since the cluster masses and the tidal 
radii are kept constant for different models, the difference in 
velocity profiles is related to the difference in cluster struc- 
ture. For example, the Aa model generally exhibits lower 
central velocity dispersion than the isotropic model. This 
means that the anisotropic model has larger r^ than the 
isotropic model. 

Near the tidal boundary, the projected velocity pro- 
files of the model Aa deviate significantly from the isotropic 
model. There are even mild bumps at around 0.5 rt for low 
mass stars in late phase (i.e., M/Mo — 0.2). Such a be- 
haviour is only seen for low mass stars in Aa models. Closer 
examination of tangential and radial velocity dispersions 
revealed that the bump is caused by the tangential com- 
ponent. The radial component monotonically decreases like 
the isotropic model. The bumps appear in almost all mod- 
els with apocenter criterion for low mass stars during late 
postcoUapse phase, but not in other models. Therefore, this 
phenomena is clearly related to the boundary conditions. 
We note here that the detailed observations by 



Drukier et al. (1998) has discovered nearly flat or 
slightly rising velocity dispersions from Ml 5 at the outer 
parts. This seems to be the indication of tangentially 
anisotropic velocity distribution in the outer parts of this 
cluster. It would be very desirable to carry out velocity mea- 
surements for other clusters. 



8 SUMMARY 

We have studied the evolution of multi-mass star clusters 
in the Galactic tidal field including the effects of velocity 
anisotropy. The radial anisotropy developed by dynamical 
relaxation tends to be suppressed by the presence of the 
tidal boundary. Except for very early epochs, high mass 
stars show more radial (larger /3) anisotropy than low mass 
stars in general. As the cluster loses a large amount of mass 
from the tidal boundary, f3 decreases, i.e., anisotropy be- 
comes more tangential for all mass species. As a result, at 
late epochs of the cluster's life, low mass stars generally have 
tangential (negative /3) anisotropy throughout the cluster, 
while high mass stars show weak radial anisotropy in the 
inner parts and tangential anisotropy in the outer parts. 
The overall degree of anisotropy depends on the half-mass 
concentration Ch of the cluster. Larger Ch clusters allow the 
development of stronger radial anisotropy in the halo. De- 
pending on the treatment of the tidal boundary condition, 
the detailed behaviour of the degree of anisotropy over ra- 
dius changes near the tidal radius. The apocenter criterion 
gives nearly tangential motion for the stars near the tidal 
radius but the energy criterion forces isotropy there. 

The density profiles of tidally limited clusters are com- 
pared with the King models. Because of complex behaviour 
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of the degree of anisotropy, finding the best-fitting King 
models becomes very difficult. The isotropic King models 
tend to give flatter profiles for low mass components com- 
pared to the Fokker-Planck results when the overall density 
profiles are matched. Introduction of radial anisotropy gives 
only worse fits, because most of the anisotropy is tangential. 

We discussed the adequacy of converting the observed 
local mass function to the global mass function assuming 
that the cluster can be approximated by isotropic multi- 
mass King models. Because of the discrepancy between King 
models and Fokker-Planck results, this process gives some- 
what erroneous results of mass function for low mass com- 
ponents. We also looked for the best place where the local 
mass function is close to the global mass function and found 
that the half-mass radius should be a reasonable place. 

The mass function inevitably changes with time because 
of the selective loss of mass. The power-law indices of the 
mass function is found to be well correlated with the frac- 
tion of mass evaporated via tidal overflow. The amount of 
fractional mass loss can be determined by observed value 
of trh and assuming the constant escape probability and 
coeval formation of all globular clusters. The observational 
data for the mass function are found to be consistent with 
the notion that the mass function is strongly influenced by 
the dynamical processes (e.g.. Richer et al. 1990; Piotto & 
Zoccali 1999), but IMF itself may have been quite different 
from cluster to cluster: some clusters may have formed with 
very steep IMF. 

The projected velocity proflles of anisotropic models 
with apocenter criterion show significant deviation from 
isotropic models or anisotropic models with energy condi- 
tion in the sense that the velocity profiles are much fiatter 
for Aa models in the outer parts during the postcoUapse 
phase. This phenomena appears to be consistent with the 
observational data for M15. 
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